Load the R packages
library(ggplot2) # graphing
library(tidyr) # dataframe manipulation
library(lubridate) # for handling date times
library(plotly) # for making interactive graphs
Read in data from EDI repository
# Package ID: knb-lter-hbr.59.10 Cataloging System:https://pasta.edirepository.org.
# Data set title: Hubbard Brook Experimental Forest: Daily Temperature Record, 1955 - present.
inUrl1 <- "https://pasta.lternet.edu/package/data/eml/knb-lter-hbr/59/10/9723086870f14b48409869f6c06d6aa8"
infile1 <- tempfile()
try(download.file(inUrl1,infile1,method="curl"))
if (is.na(file.size(infile1))) download.file(inUrl1,infile1,method="auto")
air <-read.csv(infile1,header=F ,skip=1 ,sep="," ,
col.names=c("date",
"STA","MAX","MIN","AVE","Flag"),check.names=TRUE)
unlink(infile1)
# We have a data column, but its not formatted as a date
air$DATE<-ymd(air$date) # change how R interprets Date to be a date
air$Year<-year(air$DATE)
packageID<-"Package ID: knb-lter-hbr.59.10"
Calculate annual air temp for each water year
# Use water year starting June 1st
w.year <- as.numeric(format(air$DATE, "%Y"))
june.july.sept <- as.numeric(format(air$DATE, "%m")) < 6
w.year[june.july.sept] <- w.year[june.july.sept] - 1
air$wyear<-w.year # add water year as a column to precip dataset
#subset to just HQ
hq<-air[air$STA=="HQ",]
# make sure you are only using complete years for the record
pre.obs<-as.data.frame(table(hq$STA , hq$wyear))
pre.obs[pre.obs$Freq<360, "Use"]<-"incomplete wyear" # incomplete is < 360 days
pre.obs[is.na(pre.obs$Use),"Use"]<-"complete"
# match by water year-station, only use 'complete' years
pre.obs$wys<-paste(pre.obs$Var2, pre.obs$Var1)
hq$wys<-paste(hq$wyear, hq$STA)
hq$Use<-pre.obs$Use[match(hq$wys, pre.obs$wys)]
hq.complete<-hq[hq$Use=="complete",]
# aggregate to annual average
hqa<-aggregate(list(MAX=hq.complete$MAX, MIN=hq.complete$MIN, AVE=hq.complete$AVE), by=list(wyear=hq.complete$wyear), FUN="mean")
hq_g<-gather(hqa, "type","value", 2:4)
## calculate difference between recent 6 years compared to first 6 years
magmin<-format(round(mean(tail(hqa$MIN))-mean(head(hqa$MIN)), digits=0))
magave<-format(round(mean(tail(hqa$AVE))-mean(head(hqa$AVE)), digits=0))
magmax<-format(round(mean(tail(hqa$MAX))-mean(head(hqa$MAX)), digits=0))
# order temp types
hq_g$type<-factor(hq_g$type, levels=c("MAX","AVE","MIN"))
graph air temperature at the Hubbard Brook Headquarters
## make a ggplot graph
a1<-ggplot(hq_g, aes(x=wyear, y=value, col=type))+geom_point()+ geom_line()+
ylab("Air temperature (C)")+xlab("Water year (June 1st)")+
theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
labs(col='')+
scale_color_manual(values=c(MIN="blue",AVE="black",MAX="orange"))+
geom_smooth(method="lm", se=F)+
theme(axis.text.y=element_text(size=15, colour = "black"),
axis.title.y=element_text(size=17, colour = "black"),
axis.text.x=element_text(size=15, colour = "black"))+
geom_text(aes(x=1975, y=14.5, label = paste0("Average daily high ",magmax, " degree hotter")),
size = 4.5, color="orange")+
geom_text(aes(x=1972, y=8.2, label = paste0("On average ",magave, " degrees hotter")),
size = 4.5, color='black')+
geom_text(aes(x=1975, y=3, label = paste0("Average daily low ",magmin, " degrees hotter")),
size = 4.5, color='blue')+
ggtitle("Air temperature at the Hubbard Brook Experimental Forest Headquarters")
## plotly object
p1<-ggplotly(a1)
p1
Write a html file
# this line writes the html file to create interactive graphs for the online book
htmlwidgets::saveWidget(as_widget(p1), "HQ_air_temp.html")
hubbardbrook.github.io/HQ_air_temp.html